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ABSTRACT 

We carry out adaptive mesh refinement (AMR) cosmological simulations of Milky- Way mass halos 
in order to investigate the formation of disk-like galaxies in a A-dominated Cold Dark Matter model. 
We evolve a suite of five halos to z = and find gaseous-disk formation in all; however, in agreement 
with previous SPH simulations (that did not include a subgrid feedback model), the rotation curves of 
all halos are centrally peaked due to a massive spheroidal component. Our standard model includes 
radiative cooling and star formation, but no feedback. We further investigate this angular momentum 
problem by systematically modifying various simulation parameters including: (i) spatial resolution, 
ranging from 1700 to 212 pc; (ii) an additional pressure component to ensure that the Jeans length is 
always resolved; (iii) low star formation efficiency, going down to 0.1%; (iv) fixed physical resolution 
as opposed to comoving resolution; (v) a supernova feedback model which injects thermal energy to 
the local cell; and (vi) a subgrid feedback model which suppresses cooling in the immediate vicinity 
of a star formation event. Of all of these, we find that only the last (cooling suppression) has any 
impact on the massive spheroidal component. In particular, a simulation with cooling suppression 
and feedback results in a rotation curve that, while still peaked, is considerably reduced from our 
standard runs. 

Subject headings: galaxies: formation - galaxies: evolution - methods: numerical - hydrodynamics 



1. INTRODUCTION 

In cosmologies dominated by Cold Dark Matter 
(CDM), galaxy rotation is produced by gravitational 
tidal torques arising from the hierarchical collapse of 
structure. Analytic models and N-body simulations have 
shown that this can produce enough angular mo men 
tum to ex plain the observed sizes o f disk galaxies (Fall 
fc Efstathi ou 1980; iMo et al.||1998[ ). However, compu- 
tational models including gas dynamics have struggled 
to reproduce realistic disk galaxies. Such models ini- 
tially pro duced undersized disks with low angul ar mo - 
mentum ( Navarro fc Benz|1991 Navarro & White 1994). 
Later work did genera te disks w ith approximat ely the 
correct extent ( ,Steinmetz h Navarro, 2002; Abadi et al. 
2003 1 , but these had oversized stellar spheroidal compo- 
nents and therefore unna turally large core circular veloc- 
ities (Steinmetz & Navarro 1999) far in excess of those 
found observationally. This angular momentum prob- 
lem remains one of the major shortcomings of the CDM 
paradigm. 

The origin of the angular momentum problem is not 
entirely understood, but it probably stems from the fact 
that such simulations do not achieve sufficient spatial 
and mass resolution to correctly model the appropriate 
physical processes. For example, insufficient spatial reso- 
lution leads to the spurious mixture of hot and cold com- 
ponents of the ISM, producing an artificial warm com- 

onent which is ver y efficient a t radiating away energy 

Katz||1992| |Katz et al. 1996; Ste inmetz fc Muller||1995| . 

'hus, underresolved gas in such simulations can cool very 
quickly, and in cooling it loses its pressure support and 
collapses into dense knots of material. These knots in- 
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teract with the galactic dark matter component through 
dynamical friction processes, and much of the angular 
momentum of the gas k nots is transferre d to the dark 
matter halo of a galaxy (D'onghia et al. 12006). Conse- 



quently, the knots tumble into the center ot the galaxy 
to produce a dense cusp of material in the core of the 
simulated galaxy. Typically the buildup of these mas- 
sive cores with cold gas and stars occurs rapidly, and 
even nascent simulated galaxies exhibit evidence of these 
cusps as early as z ~ 4-5. 

Energetic feedback from star formation events has the 
ability, in theory, to alleviate this problem either by heat- 
ing and "puffing" up the collapsing knots so they are 
more easily disrupted b efore losing their angular momen- 
tum (Weil et al. 1998), and/or by preferentially eject- 
ing low angular niomentum gas (e.g., ,Brook et al._2010J . 
However, existing simulations are unable to resolve the 
detailed structures of star-forming regions. Individual 
star-forming and stellar feedback events occur on parsec 
scales. Stellar feedback processes can be resolved in sim- 
ulations confi ned to local regions of the ISM (e.g., Joung 
fc Low|[2006 ) but not in cosmological simulations, which 



need to accurately co-evolve a galaxy's environment on 
scales of > 10 Mpc. 

Therefore, cosmological simulations simplify and pa- 
rameterize star formation and stellar feedback on scales 
more easily met by current computational resources (i.e. 
scales of 10^ — 10"^ pc). A variety of techniques have been 
suggested in order to achieve this puffiness within the 
confines of low-resolution models. In the most primitive 
prescription, stellar feedback is simply the return of en- 
ergy from newly created star particles to their surround- 
ing gas, usually in the form of thermal energy. Star "par- 
ticles" in these simulations typically represent clusters of 



2 



stars of mass 10^ — 10^ M0, so the thermal feedback 
is justified as the energy output from the most massive 
stars in the cluster becoming type II supernovae shortly 
after the creation of the star particle. In nature, this 
supernovae heating produces a small component of very 
hot gas surrounding the stellar population, so hot that 
its cooling time is very long. Unfortunately at the res- 
olution currently achievable in cosmological simulations, 
this primitive thermal prescription for feedback deposits 
the same SN energy into a much larger reservoir of gas, 
which does not reach the same high temperature as it 
should. This now-warm component of gas can easily ra- 
diate the excess energy away, cool further and proceed 
with star-formation r unaway thus defeating the p urpose 
of the feedback (e.g., |Stei nmetz fc Navarro||1999 ). 

Building from these failures, a number of research 
groups artificially turn off cooling in a gas parcel for a 
period of time {t ^ 10'' yr) afte r a cluster o f stars has 
formed out of it (e.g., Gerritsen|l 997; Thacke r fc Couch- 



stiffer in order to provide an additiona l source of pressure 



man 



2006 



2000^ 'So mmer-Larsen et al.||2003( iStinson et al.' 
Governato et al.||2007[ [Agertz et al.||20101 IColi: 



eTatlTSOlOj iPiontek fc Steinmetz||2011[ |Guedes et al. 



2011 1. This method is justified as an application of the 



Ke dov- Taylor blast wav e solution for a Type II SN ( Tay 



lor 1950 Sedov 1959), which blows out any cold me- 



dia from the immediate environment of a star formation 
event. Using this prescription, any gas in a galaxy which 
starts to collapse into knots will reach the star forma- 
tion criteria, form a star, and then heat up without any 
allowed cooling, thus preventing further collapse. Not 
surprisingly, these research groups have found some suc- 
cess with this method, yielding simulated galaxies with 
reduced inner rotation curves due to less massive bulge 
components; however, gas parcel masses and sizes in cos- 
mological simulations of this sort are typically to o larg e 
for the Sedov- Taylor solution to apply (see Section 2.1.7). 



Thus despite the successes of the cooling suppression 
feedback model, the community continues to search for 
other more physically-motivated solutions. 

Another subgrid model for feedback (i.e. on scales 
smaller than the true resolution of the simulation) is to 
inject kinetic energy directly into the gas; this can allevi- 
ate the problem of thermal energ y being radiated away. 



For example, some studie s (e.g. 'Springel & Hernquist 
Scannapieco et al.| [2006 ; Oppcnhcinicr & Dave 



2003 



2008p using Smoothed Particle Hydrodynamics (SPHj 
give some of the SN energy to individual gas particles 
in the form of momentum. This method can result in 
significant mass outfiows (by design), but at the cost of 
decoupling wind particles from hydrodynamic interaction 
for a period of time. An alternate approach, to keep wind 
particles coupled to the disk gas was explored by Schaye[ 
& Dalla Vecchia (2008). Both approaches help but, Ey 



themselves, do not appear to generate realistic rotation 
curves. 



In addition, Truelove et al. (1997) showed that insuf- 
ficient resolution in a simulation can lead to artificial 
fragmentation of the gas, perhaps resulting in a further 
overproduction of stars. One way to prevent artificial 
fragmentation is to add additional (numerical) pressure 
in high-density, low-temperature r egions to ensure that 



the Jeans length is always reso lved (Machacek et al. 2001 
Robertson & Kravtsov[ 2008[). This can be achieved by 



[R^ 

modifying the equation o± state (EOS) itself, making it 



to gas m denser regions ( bcnaye 6z JJ alla veccnia /UUo 
Ceverino fc Klypin[p009| [Agertz et^l.||2010i ). A poly 
tropic EOS (P a p^ ) with P = 4/3 will keep the ratio 



of Jeans length to resolution length constant (assuming 
Lagrangian resolution such that the resolution length de- 
creases as p~^l^ - for fixed resolution P = 2 is required), 
but even stiffer rela tions have been used. For example. 



Agertz et al. (20101 ran simulations with such an equa- 
tion ot state, where in low-density regions it behaved as 
an ideal gas, but in high-density (star forming) regions 
it followed a polytropic equation of state with P = 2. 

In this paper, we undertake an investigation of galaxy 
formation using an Adaptive Mesh Refinement (AMR) 
hydrodynamics code. The majority of work in this field 
has used SPH codes, and so this allows us to investi- 
gate the problem from a new angle. Although there has 
been some work with AMR codes (Joung et al. 2009a[ 
Ceverino fc Klypin||2009l [Agertz et ai.||20i0[ |Colin et al' 
2010[ ), there has not been a clear demonstration that an 
equivalent AMR calculation (i.e. one without a subgrid 
feedback model) actually does reproduce the classic SPH 
result. 

We begin by simulating a set of five halos without any 
feedback or subgrid model (except a minimum pressure 
support to prevent artificial fragmentation). We find, 
in agreement with SPH codes that a large, concentrated 
bulge is produced, resulting in a rotation curve that rises 
to ~ 500 km/s at r ~ 1 kpc. We then vary a number 
of numerical and physical parameters in order to under- 
stand how sensitive the result is to our a choice of pa- 
rameters. 

The paper is organized as follows. Section [2[ describes 
the details of our hydrodynamics code, our initial con- 
ditions and the relevant parameters for this study. In 
Section [3[ we present the results of our simulations in- 
cluding the five canonical runs, our resolution study, and 
our modified runs. Scction[4|is a discussion of our results 
and their implications. Finally, Section[5]summarizes our 
conclusions and makes predictions for future solutions to 
the angular momentum problem. 

2. METHODOLOGY 
2.1. Code 

Our simulations were performed using iSnzcj^ an Eule- 
rian, three-dimensional, grid-based hydrodynamics code 
that employs adaptive mesh refinement in order to 
achieve targ eted regions of h igh resolu tion in a cosmolog - 
ical volume ( [Bryan fc Norman 1997; O'Shea et al.|2004p . 
Gas is discretized on the grid, but dark matter and stars 
are treated as particles. The ZEUS hydrodynamics code 
( Stone fc Norman[[1992 1 is used to evolve the gas on the 
grid. Unzo includes gas, self-gravity, a non-equilibrium 
model for H and He ionization a nd cooling, a metagalac- 
tic ultraviolet background (Haardt fc Madau 1996| ) , and 
equilibrium cooling due to metals (although tor the runs 
described in this paper, we do not include metal cooling). 



2.1.1. Star Formation 
Star formation is modeled using a simple criteria based 



on 



Cen & Ostriker ( 1992 ). A grid cell will produce a star 



enzo-pro ject . org 
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if: (i) the overdensity in that cell exceeds a given value 
((Ssf), (ii) the mass of gas in the cell exceeds the local 
Jeans mass, (iii) there is locally convergent flow (i.e. the 
velocity divergence is negative) and (iv) the cooling time 
for the gas to collapse is less than the dynamical time 
in that cell (or the temperature is near the minimum 
allowed, around 10'* K). If a grid cell meets all the previ- 
ous criteria then gas is converted into a "star particle", 
as calculated using 



At 

eSFT Pgas 

^dyn 



(1) 



where esp is the star formation efficiency (more properly 
the efficiency per dynamical time). At is the size of the 
time step, idyn = (37r/32Gp)^/^ is the dynamical time 
and pgas is the gas density. If the resulting star particle 
would have a mass above a minimum mass, uiiiu it 
is formed immediately. This mass is chosen so that a 
large number of small star particles will not slow down 
the simulation. However, if the star particle would have 
a smaller mass, the probability that it will form is equal 
to the ratio of the mass of the projected star particle 
to M*^„iin. If the star particle is then formed, its mass 
is the minimum of M^.niin and 90% of the mass in the 



gas cell (see also Tasker & Bryan 20061. For all of the 
simulations in this study, we used a value of M^^min = 10^ 
M0. Our default value for esF = 10~^, although this was 
varied in some runs (see Tablc[l]). We adopt ^sf = 1000, 
corresponding approximately to a density threshold of 
0.1 cm~^ at z = 3; this value simply ensures that star 
formation is limited to strongly non-linear regions. 

2.1.2. Resolution 

Because Enzo is an AMR code, it can dynamically re- 
fine the spatial resolution of the grid when certain criteria 
are met in a particular region of the simulated volume. 
We set the grid-refinement criteria to increase refinement 
whenever the dark matter mass in a cell in larger than 
four times the dark matter particle mass, with an equiv- 
alent criterion for the baryonic mass. When refined, the 
cell resolution is increased by a factor of 2. The place- 
ment of these refinement regions is recalculated regu- 
larly throughout the simulation, to assure that moving 
or emerging regions of interest are always well-resolved. 
For the canonical runs, we cap this refinement when our 
cell sizes reach 9 levels of refinement or 425 comoving 
parsecs. In a ddition, we conduct a resolution study in 
Section 3.5.1 where we vary the resolution from 7 to 10 
levels of refinement (212 - 1700 comoving parsecs). 

2.1.3. Minimum Pressure Support 

To prevent artificial Jeans fragmentation, for most runs 
w e implement the minim um pressure support described 



Machacek et al. 



Jeans length to ce 



(2001|) such that the ratio of the the 
Lj/Ax, is at least 



'1 size, J = Lj/A^, is at least 8. We 
do this by adding an additional artificial pressure to the 
most highly refined grid cells such that this ratio is always 
maintained. The addition of this pressure is intended to 
prevent gas clouds from collapsing below the resolution 
scale of the simulation, which could cause spurious nu- 
merical effects, such as artificial fragmentation. 
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Fig. 1. — We plot gas surface density versus star formation rate 
surface density for our five canonical runs. We overplot the power 
law SsFR = 2.5 X fO~* ■ ^gas (in these units) representing the 
Kennicutt-Schmidt local relation bet ween star formation rate and 
gas surface density l |Kennicutt|1989l l. 

2.1.4. Star Formation Efficiency 

For our canonical runs, we set the star formation effi- 
ciency per dynamical time esF to 1%, a parameter value 
that previous work found to approximately reproduce the 
Kenn icutt-Schmidt relation for our chosen star formation 
law (Tasker fc Bryan|2006[ ) . In Figure [l] we demonstrate 
that our canonical simulations roughly agree with the 
Kennicutt-Schmidt relation by plotting the gas surface 
density versus the recent star formation surface density 
(agCstar ^ 5 Myr) in concentric annular cylinders of 0.25 
kpc width and scale height 5 kpc each aligned with their 
respective gas disks. Overplot is a power-law with in- 
dex 1.4, the accepted form of the Kennicutt-Schmidt law 



(Kennicutt 1989). Our canonical value for the param- 
eter (the star formation efficiency per free fall time) is 
esF = 10~^, which is in line with ( but slightly lower 
than) typically suggested values (e.g. Krumholz fc Tan| 
2007 1. Other galaxy for mation simul a tions adopt ditter- 
ent valu es, for example Abadi et al. (200 3| use a value 
of 3.3% , Governato et aTl2U07) use 57o, iStinson et al. 



( |2006| use 5 %, all tor SFH simuTations. Turning to AMK 
simulations, Agertz et al . (2010) used 1% for their most 
successful runs, I'l'eyssier et al.| |2010l also adopted 1%, 



while Gnedin fc Kravtsovl2Ull ) adopted 0.5%, arguing 
that it was a better fit to rhe observations of [Bigiel et al.| 
( 2008 1 . It is likely that the choice of star formation ef- 



ficiency required depends both on the resolution of the 
simulation, as well as feedback prescription. In this pa- 
per, we explore variations in the efficiency, decreasing 
this parameter to explore the impact of slowing the con- 
version of gas into stars. 

2.1.5. Fixed Physical Resolution 

Most of the previous work on this topic was done us- 
ing hydrodynamics codes in physical or proper units as 
opposed to comoving units. Our code uses comoving co- 
ordinates, which means that a fixed maximum resolution 
level therefore corresponds to a fixed comoving resolution 
at the finest refinement level (cold, dense gas is almost 
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always refined at the maximum level in these simula- 
tions) . Thus our models resolve more details in the early 
universe (by a factor of 1 + z) relative to simulations em- 
ploying the equivalent resolution in fixed physical units. 
In an effort to explore the impact of using fixed physical 
rather than comoving resolution, we modified the code 
to be able to run with a fixed physical resolution, by ad- 
justing the maximum refinement level with redshift. We 
note that because of the discrete, factor of two changes in 
grid resolution, this only keeps the maximum refinement 
to a fixed physical resolution within a factor of two. We 
do not change the refinement criteria, ther efore our mass 
resolu tion remains the same (see Section 2.2). Section 
13.5.41 describes the results of these runs. 

2.1.6. Supernovae Feedback 

Most simulations we discuss in this paper do not in- 
clude explicit feedback; however, we do carry out a few 
runs which included prompt, energetic feedback from 
Type II SNe. We apply a simple prescription for stellar 
feedback to mimic the effects of type II supernovae. Be- 
cause individual star- format ion events produce star par- 
ticles of Mgtar 10^~^ Mq, we spread feedback over an 
extended period parameterized by 



dM, 



form 



dt 



Mn 



exp 



(2) 



where Mq and tg are the initial star particle mass and 
star particle creation time (see also Ccn fc Ostriker^l992j ), 
and T is the maximum of either the dynamical time ot the 
gas out of which the star particle formed, or 10 Myr (to 
prevent unrealistically short dynamical times). We take 
a fraction of the rest mass energy epB of the forming 
stars as the available feedback energy. This parameter 
can be computed assuming an initial mass function and 
a minimum initial stellar mass for producing a Type II 
SN. 



As described in Tasker & Bryan (2006), we simply add 
the thermal energy to the local grid cell. As the sur- 
rounding gas heats up, it increases the Jeans length and 
theoretically quenches star formation for a time. As was 
noted in Section [T] if the resolution of the simulation is 
insufficient to clearly differentiate between a cold, neu- 
tral interstellar medium and a hot, ionized interstellar 
medium, then the pumping of this energy into the cold 
star-forming gas results in an unrealistic warm compo- 
nent (i.e. the mixing of hot and cold phases). This warm 
component sits near the peak of the cooling curve of the 
gas, and thus effectively radiates away its energy very 
quickly. 

For the run which does include feedback, we take a 
value of epB = 3 x 10~^, which corresponds to one 10^^ 
erg SN for every 180 Mq of stars formed. The value of 
this parameter is uncertain, as it depends on both the 
initial mass function, as well as assumptions about how 
the energy is radiated away immedi ately. Other values 
are used in the literature, for example Abadi et al. ( 2003 1 
employed a feedback prescription which corresponds to, 
in our definition, a value of epB = 5.6 x 10~^. They, 
like others, found that this energy was quickly radiated 
away, and so the value ha d little impact . As we dis- 



cuss in more detail below, Stinson et al. (2006) added 
a cooling suppression model, and argued tor a value 



of epB = 4.3 X 10"'', although [Covernato et al.| ( |2007 l 
adopted epB = 2.6 x 10"^, using a very similar model. 
Turni ng from SPH to AMR simulations, we note that 
Agertz et al.|(|2010) used efb = 5.6 x 10"^ for their stan- 
dard runs (also with a co olin g suppression model). Nei- 
ther [Teyssier et'al] ( |2010l ) nor |Gnedin fc Kravtsov| ( |2011 1 
appear to hav e "used ther mal feedback in their simula- 
tions. Finally, |Cen| ( |201l| used epB = 10 x 10"^ (larger 
than the usual etticiency because of a postulated contri- 
bution from prompt Type la SN); this work did manage 
to drive winds without a cooling suppression model, al- 
though they added the energy to the nearest 27 cells 
weighted inversely by density (and also had somewhat 
worse mass and spatial resolution than used here). In 
summary, we see that our chosen value is within the range 
used by other researchers. 

2.1.7. Cooling Suppression Models 

One way to prevent this energy from being quickly 
lost is to turn off radiative cooling in the region immedi- 
ately surro unding newly for med stars. This was first at- 
tempted by Gerritsen ( 1997 ) , but has since been explored 
by a range ot simulations (e.g., Thacker fc Couchman 



2000| |Sommer-Larsen et al.||2U03| |Stlnson ei al.limi 



Governato et al.||2007[ [Agertz et al.||2010| |CoKn et al.| 



2010| [Piontek fc Steinmetz||2011l [Guedes et al.||201ip . 

I'he idea is to use the Sedov-'l'aylor solution for a blast 



ayli 

wave to model the subgrid shock physics that the code 
cannot resolve. It might seem most straightforward to 
use the length and time-scales of the energy-conserving 
phase (i.e. Sedov phase) to control where and when to 
turn cooling off, but these turn out to be too small and 
too short (a few 10'' years) to make a sig nificant differ 



ence, as acknowledged by previous studies ( Stinson et al 



2006 1 . Instead, the method employs the radius and time 



at the end of the momentum-conserving phase under the 
assumption that during this phase, much of the energy is 
conserved in kinetic motion (and hot, diffuse gas), which 
the code cannot model, and would otherwise dissipate. 
The larger length and timescales of the snowplow phase 
used are further justified as the combined forces of mul- 
tiple supernovae in the star particle; however, although 
these supernovae may interact in a complex way, it is 
unlikely their effects will simply add constructively. 

This prescription for cooling suppression feedback has 
generally been used in SPH codes, but recentl y some 
AMR codes have also adopted this technique (Agertz] 
jet al. 2010; Colin et al. 2010). The common prescription 
IS to suppress cooling for a period of time (30 — 50 Myr) 
in the gas immediately around a star- formation event, 
regardless of where the star goes afterwards. In our im- 
plementation, we suppress cooling of the gas in the single 
cell in which the star particle resides. This is done for 50 
Myr after the star particle is first created. Since both of 
these length and time scales correspond closely to those 
over which energy from the star particle is injected in 
the simulation (the feedback follows Eq. [2] above) , this 
acts to suppress cooling in newly heated gas. Given our 
chosen cell size (425 comoving pc in most runs), these 
length and timescales are similar to the region and du 



ration of influence adopted by other researchers ( Stinson 
eraL][2006l [Coh'n et al.||2010| . 



2.2. Initial Conditions 
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''Halos H26SPM and D9H26SPM are the same simulation. 

''Halos H26SPMR and H26SPMFCR have a fixed physical resolution throughout the simulation, while all other halos have a fixed comoving 
resolution. 



TABLE 1 

We present the bulk properties and simulation parameters for each halo and simulation. The first five halos are the 
canonical simulations, each following a different halo in the same overall volume. additional simulations were created 
using the initial conditions for halo h26spm with modifications to the encoded baryonic and star formation physics. 



For tliis work, we use t he WMAP 5-year resul ts as our 
cosmological parameters (Komatsu et al. 2009), in par- 
ticular we use: fla = 0.258, Ha = U.742, iibaryon = 0.044, 
<Ts — 0.796, Ho = 71.9 km s^^. We generate our initial 
conditions using inits, a program included in the Enzo 
suite. Inits sets up a 128^ particle-mesh grid, and modi- 
fies the velocity and position of the dark matter particles 
in each grid cell as specified by linear perturbations with 
the required power spectrum at z = 99. The initial con- 
ditions are generated for a cubic volume oi L ~ 20h^^ 
comoving Mpc on a side with periodic boundary condi- 
tions. First, the simulation is populated with only dark 
matter particles at low-resolution (Mdm — 3.2 x 10* Mq) 
and run to z = 0, where candidate Milky- Way-like halos 
are identified. Halos are chosen based on their final mass 
and accretion history, preferentially selecting halos with 
final masses M2Q0 ^10^^ JVI©, and those which have not 
undergone major mergers after z 2 (M200 is the mass 
enclosed within a radius corresponding to a mean den- 
sity of 200 times the critical density). Five such halos 
are identified, ranging in mass from M200 = 4.8 x 10^^ 
to 1.2 X 10^2 

For each halo, the component dark matter particles 
are traced back to their positions in the initial condi- 
tions at z = 99. This Lagrangian volume is further re- 
fined with two additional levels of refinement. It is here 
that the initial conditions are regenerated, and in these 
nested boxes, we additionally refine the dark matter par- 
ticle masses by a factor of 8 for each region. The result- 
ing high-resolution dark matter particle mass within the 
vicinity of each halo is Mdm = 4.9 x 10^ M0. 

A series of new simulations are performed using these 
new initial conditions; each one focuses on a different 
halo. Baryons are included in these runs. The high- 
resolution regions are further refined dynamically with 
adaptively-placed grids, using the refinement scheme de- 
scribed earlier. 



2.3. Description of Simulations 

We conducted five canonical simulations using iden- 
tical simulation parameters and the initial conditions 
described above. These simulations are referred to as: 
H26SPM, H30SPM, H37SPM, H47SPM and H54SPM. 
Additionally, several different runs were performed on 
the initial conditions for halo 26 (H26SPM), which sys- 
tematically varied simulation and physical parameters 
to investigate the effects of each on galactic evolution. 
The parameters toggled (and their respective simula- 
tions) include: (i) excluding minimum pressure support 
H26S]; (ii) changing the maximum spatial resolution 
D7H26SPM, D8H26SPM, D10H26SPM]; (iii) using a 
constant physical resolution instead of a constant comov- 
ing resolution [H26SPMR, H26SP1VIFCR]; (iv) including 
thermal feedback [H26SPMF, H26SPMFCR]; (v) lower- 
ing the star- format ion efficiency [H26SPML]; and (vi) 
suppressing cooling in star forming regions [H26SPMC, 
H26SPMFC, H26SPMFCR]. The details of the various 
simulations and the resulting galaxies are shown in Ta- 
bled! 

3. RESULTS 

In this section, we present the results of our galaxy 
formation simulations, first describing the five canonical 
runs, which all contain identical physical prescriptions 
but track different galactic halos. Then, we explore vari- 
ations in resolution as well as the numerical parameters 
we use to describe the gas and star formation. 

3.1. Mass History 

The mass accretion history for a galaxy including the 
different modes of its accretion is thought to play a cru- 
cial role in deter mining its final dynamical state (e.g. 
Keres et al.||2005 ). In order to analyze the simulation in 
high time resolution, we record outputs from the simu- 
lation ever y 10 Myr. For each ou tput we run the HOP 



algorithm (Eisenstein & Hut 1998, ) on the dark matter 
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Fig. 2. — We plot the mass of our five canonical halos as a function of time. Different components of each halo are labeled differently: 
dark matter (green dashed line), stars (red dot-dashed line), gas (blue dotted line), and total mass (black solid line). 



particles in order to identify halos. Given the particles in 
each of our five halos at z = 0, we identify and track these 
halos back to early times. Each halo is tracked backwards 
in time by identifying the local progenitor which shares 
the largest number of tightly-bound dark matter parti- 
cles. The resulting mass-accretion history for each halo 
is shown in Figure |2j The halo masses are computed in- 
side of r2oo, the radius within which the mean density is 
200 times the critical density of the universe at that red- 
shift. At each time, we determine the center of the halo 
using an iterated center-of-mass technique, which starts 
with the center of mass within r2oo j and then successively 
recomputes the center of mass in smaller spherical vol- 
umes, decreasing the radius by 5% on each iteration and 
using the center of mass of the previous volume. This is 
necessary in order to make an accurate determination of 
the halo center (we found that simply choosing either the 
densest point or the center of mass within r2oo did not 
produce a good estimate of the center in many cases). 
All masses calculated are masses contained within r2oo, 
and are shown in Table [l] 

In Figure [2) the mass histories of dark matter (green) , 
gas (blue), stars (red) and total mass (black) are shown. 
The halos are arranged in decreasing z=0 total mass. 
All of the halos lack any major mergers over the last 
10 Gyr, providing them with quiescent growth, in order 
to maximize the chance of producing disk systems. The 
lowest-mass system, halo 54, undergoes mergers at z ^ 
3.5 and z ~ 2.5 causing discrete jumps in the mass of all 
its components at those times. 

3.2. Star Formation History 

Much of the angular momentum problem stems from 
an overproduction of stars and a buildup of the oversized 



stellar bulge, revealed by the galactic star formation rate 
history. If one can reduce star formation early in a halo's 
evolution it will moderate the amount of material in the 
inner region of the galaxy. While it is true that the bulk 
of the modeled bulge mass is stellar in nature, it remains 
unclear as to where these stars were created. Were they 
formed in clumps in the disk before plummeting to the 
center of the system, or did dense knots of gas spiral into 
the center of the galaxy where they ultimately collected 
and formed stars ? This question will be examined more 
closely in Section f3.5.3[ 

In Figure [31 we show the star formation histories for 
our halos. These are computed by selecting the particles 
within r2oo at z = and using the stellar age of each par- 
ticle to compute the implied star formation rate. This 
means that the rate shown includes star formation in all 
progenitor halos (not just the most massive progenitor 
shown in Figure 2]) . The star formation history shows an 
early burst, as dense, cold gas is rapidly converted into 
stars, followed by a decline to a steady continuous level of 
SFR ~ 3 — 10 M0 yr~^. This low- level of star formation 
reflects both the decreasing gas suppl y and the dec rcas- 
ing amount of cold gas accretion (e.g. Keres et al.|2005) . 
Its value is slightly higher than but roughly consistent 
with observed levels in the M ilky Way of ~ 1 M0 yr~^ 
(Robitaille & Whitney 2010). Interestingly, our lowest 
mass halo, H54SFM, undergoes a much less pronounced 
early burst of star-formation, never forming more than 
15 MQyr~^, which may be a reflection of a more extended 
merging period in the first few gigayears of evolution, as 
evidenced by its accretion history in Figure [2j 

3.3. Disk Images 
3.3.1. Gas 
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Fig. 3. — We plot the star formation histories for the five canonical halos. The bulk of the star formation occurs in the first 2-3 gigayears 
of evolution for each halo before settling into a low level of continuous star formation lasting to present. 



Each of the five halos pro duce a gaseous di sk at z = 0. 
We use the analysis suite yt ([Turk et al.|2011|) in order to 
visualize the gaseous and stellar components of these ha- 
los, as shown in Figure l4] We determine the disk normal 
vector by computing the net angular momentum of all 
cold, dense gas (defined as T < 2 x 10* K and p > 10^"* 
MoMpc"'^) within 0.2 r2oo- We then generate two im- 
ages for each disk, one side-on and one faceon but each 
with the same scale of 25 kpc on a side. These two images 
are generated in slightly different ways. For faceon pro- 
jections, we simply show the gas surface density in units 
of particles per square centimeter. For the edgeon image, 
we carry out a volumetric rendering to show transparent 
isodensity contours (in detail, we use the gas density to 
assign a transfer function at each point that consists of a 
set of narrow Gaussians, each separated by a factor of 10 
in density, and each given a different color). This allows 
us to show both the disk but also bring out structure in 
the halo. 

These images display gas disks with radii of a few lO's 
of kpc, slightly smaller than, but approximately typi- 
cal of present-day late-type systems. These cold disks 
are present from early times and are rotationally sup- 
ported - we will examine their rotation al v elocities in 
more detail in the next section (Section 3.4). None of 



the disks display significant spiral structure at this par- 
ticular timestep, although all of them seem to possess it 
at some point in their evolution. Interestingly, the most 
massive of the halos, halo 26 does not display as large 



of a disk as halo 47 despite having nearly two times as 
much total mass; however Table [T] reveals that the latter 
has significantly more cold gas, which is what is plotted 
here. The details of satellites and gas accretion tend to 
dominate the observed behavior of the disk at any point 
in time, and halo 26 recently accreted some small halos. 

3.3.2. Stars 

In the bottom part of Figure |4j we render the stellar 
particle distribution for each halo. Using the same scale 
and camera angles as we did for the gaseous components 
above, we generate stellar surface density plots. We rep- 
resent each star particle as a gaussian with a sigma of 200 
pc, then step through the volume depthwise, coadding a 
random sampling of 10% of all star particles. The color 
of a star is determined by its age, where we use a contin- 
uous color function of blue through red to represent the 
log of the age of the star as shown in the colorbar. This 
represents stars just 10 Myrs old as pure blue and stars 
13 Gyrs old as pure red. Just as in the gas surface den- 
sity, the intensity of our stellar renderings is logarithmic 
with the bulk of the material residing in the inner core 
of each galaxy. 

All of the stellar halos lack the disk structure we might 
expect from disk galaxies. Instead, they are completely 
dominated by a bulge component. This happens because 
star formation occurs most efficiency in the small, dense 
clumps that merge and lose angular momentum. With- 
out feedback, the dense gas clumps efficiently deposit gas 



Fig. 4. — We render faceon and edgeon views for each of our five canonical fialos at z = 0. The top two rows show the gas component 
of each galaxy, whereas the bottom two rows display each galaxy's stellar component. Each postage stamp has a width of 25 kpc. For the 
gas, edgeon views show a volume rendering with isocontours of the gas density, whereas faceon views are column density bricks. For the 
stars, we use the column density of stars in the faceon and edgeon views respectively. The color of a star particle is representative of its 
age where blue represents young stars and red represents old stars according to the color bars on the right. 



and stars in the center of the halo, leaving only a gas-poor 
disk which does not efficiently form stars. In addition, 
we do not see significant age gradients for most of the 
halos. We note that halo 47 does display the youngest, 
bluest overall stellar population (particularly relative to 
halo 30's aging stars), consistent with the star formation 
histories of Figure [Sj 

3.3.3. Larger gaseous environment 

In addition, we generate large-scale volumetric render- 
ings of halo 26 in order to show its extended halo and 
immediate environment. Figure [5] is generated in the 
same way as the edgeon renderings described in Section 
3.3.1 It shows H26SPM over a region with 250 kpc (ap- 
proximately the virial radius) and 2.5 Mpc on a side re- 
spectively. These large-scale volume renderings (of the 
side-on disk) show that there is also a gaseous halo, ex- 
tending out to at least r2oo ^ 200 kpc. This hot halo 
gas is approximately spherical at z = 0, but contains a 
significant amount of substructure due to ongoing infall 
and asymmetric accretion. We do not see gas clumps 
coohng a nd condensing out o f the smooth, hot g as halo 



(see also Binney et al. 2009 Joung et al. 20111. The 



larger-scale image shows that the halo is embedded in a 
set of filaments, along which gas accrete and is typical of 



B 



Fig. 5. — We display edgeon volume renderings of Halo 26 at z = 
0. Like the edgeon views of the galaxies in Figure |4] these images 
show structure through the use of isodensity contours; however, 
unlike Figure |4] these images show the larger environment around 
the galaxy witn 250 kpc and 2.5 Mpc on a side respectively. 

the other halos. 

3.4. Rotation Curves 

In this section, we focus on rotation curves at z = 0, as 
these are both directly comparable to observations and 
also immediately show the mass distribution of the halos. 
In Figure [6j we plot the the rotation curves for each of 
our five systems. The curves in each graph represent 
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Fig. 6. — We plot the rotation curves for the five simulated halos. We represent the circular velocity of the gas due to the dark-matter 
component (green dashed line), the stellar component (red dot-dashed line), the gas component (blue dotted line), and the total of all halo 
components (black solid line). The stellar component seems to be driving the cusped rotation curve in the cores of each galaxy, whereas 
the gas and dark matter profiles appear consistent with observational expectations. In each system, the total-mass rotation curve is highly 
cuspy and unlike any observed galaxy, confirming that our simulations reproduce the angular momentum problem. 



the equivalent circular velocity for each mass component: 
dark matter, gas and stars: 



GAf(< 



(3) 



This figure shows that the rotation curves for each halo 
are highly peaked in their inner 5 kpc, primarily due to an 
overabundance of stars in their cores. This is clearly in- 
consistent with the nearly flat rotation curves observed in 



disk systems (e.g., Courteau][l997 ). The gas contributes 
negligibly, while the dark matter curve is steeper than ex- 
pected for an NFW-profile because of contraction driven 
by the deep potential well of the stellar component, and 
remains nearly flat to the core, but is secondary to the 
stellar distribution. 

3.5. Modifications to the Canonical Runs 

In addition to our five canonical runs presented above, 
we conducted a series of additional simulations in which 
we systematically varied numerical parameters involv- 
ing resolution, star formation, supernovae feedback, and 
gas physics. We used the initial conditions from halo 
H26SPM for each of these simulations, so that we could 
directly compare these results against one of our canon- 
ical models. In the following subsections, we present the 
results of these various modified runs and for each sim- 
ulation, examine its star-formation history as well as its 
rotation curve. The plots for all of these runs are pre- 
sented side-by-side in Figures [TOl [TT] [12] & [13| The bulk 
characteristics of these halos are presenteoin Table [T] 
We note that some runs of this halo have a slightly higher 
virial mass at z = due to the presence of a fairly large 



satellite (10% of the main halo mass), which sits very 
close to the virial radius - slight shifts in its position 
in the various runs can lead to its inclusion in the total 



3.5.1. Resolution study 

Since spatial resolution has been raised as an impor- 
tant issue affecti ng the angular momentum of the gas 
(e.g.|Ma,yei;[2004|, we conducted three different modified 
runs, each witTia factor of 2 change in the spatial reso- 
lution (done by systematically changing the maximum 
allowed refinement level). Runs D7H26SP, D8H26SP 
& D10H26SP follow the initial conditions of H26SPM 
(which itself has a resolution of 425 comoving pc) with 
1700, 850, and 212 comoving parsecs spatial resolution 
respectively. There is no change to the mass resolution 
in these modified runs. 

Figures [7] and |8] show the star formation histories and 
z — rotation curves for these runs. Up to the range we 
can probe, resolution alone doesn't appear to play a sig- 
nificant role in determining the overall mass distribution 
in a halo or the conversion of gas into stars. There is an 
indication from Figure [7] that lower resolution lowers the 
initial burst (because of the decreased central gas den- 
sities) and therefore shifts the peak of star formation to 
later times; however, the effect is small. Increasing the 
spatial resolution actually allows material to condense to 
an even smaller volume in the cores of galaxies, which 
increases the inner cusp in the rotation curves of these 
systems, as seen in Figure [8j Thus, resolution by itself 
cannot solve the angular momentum problem. Perhaps 
increased resolution when coupled with a more sophis- 
ticated star-formation or feedback prescription will pro- 
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Fig. 7. — We plot the star formation histories for our resolution study. These runs all use identical initial conditions, those of canonical 
run H26SPM, which is also presented here as D9H26SPM. The only difference in each run is the maximum refinement level, that is, the 
maximum level of spatial resolution achieved which ranges from 212 comoving parsecs to 1.7 comoving kiloparsecs. 
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Fig. 8. — We plot the rotation curves for the resolution study. As in Figure |6] we represent the circular velocity of the gas due to the 
dark-matter component (green dashed line), the stellar component (red dot-dashed line), the gas component (blue dotted line), and the 
total of all halo components (black solid line). There is a clear trend in these runs that as one increases the resolution of the simulation, 
one allows a denser cusp of material in the core (primarily in the form of stars). In turn, the core circular velocity is increasingly driven 
upward. 
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duce more realistic galaxies when one begins to resolve 
star-forming regions on parsec scales. 

3.5.2. Minimum Pressure Support 

These simulations cannot resolve star-forming events 
on parsec scales, so our canonical runs employed an arti- 
ficial minimum pressure to prevent Jeans fragmentation 
on smaller scales than we could resol ve. By includ i ng thi s 



minimum pressure described in Machacek et al. (2001), 
we assured that the Jeans length is always retined by 
at least 8 cells, and therefore the Truelove criteria was 



met (see alsofTruelovc et al. 11997 


Robertson & Kravtsov 


2008 Ceverino et al.||2010). It has been suggested that 



artihcial disk fragmentation leads to large gas clumps 
that lose angular momentum via dynamical friction, and 
hence result in angular momentum loss. For modified 
run H26S, we turned off this minimum pressure support. 
The results of that run are presented in Figure [lO] and 



11] The results are nearly identical for both measures, 
indicating that the minimum pressure support did not 
have a significant effect on the star formation rate or the 
distribution of gas and stars in our simulated system. 

3.5.3. Star Formation Ejficiency 

The star formation model we adopt in this work is 
based on the Kennicutt-Schmidt relation, but has an ef- 
ficiency parameter that is not well-constraine d. It has 
recently been suggested by Agertz et al. (j2010') that star 
formation efficiency is a key parameter controlling the 
distribution of stars in the disk, and hence the rotation 
curve. Our efficiency is already fairly low; however in 
order to examine this suggestion in more detail, and also 
to probe the impact of decreasing the efficiency in gen- 
eral, we adopt esF = 10"^ in run H26SPML. The star 
formation history, in Figure |10[ shows a significant delay 
to later time, as would be expected, although the net 
amount of stars produced is nearly identical. However, 
this does not translate into a flatter rotation curve, as 
can be seen in Figure [TT] 

To understand this result a little more, we look at the 
rotation curve for this run, at early times. Figure [9] shows 
the rotation curves aX z — 4.5 of this low star formation 
efficiency run H26SPML compared against our canonical 
run H26SPM. During this epoch, much of the halo has 
formed, but star formation has not yet converted most of 
the gas into stars. This plot demonstrates that gas dom- 
inates the rotation curve at early times, and yet despite 
this, the gas clumps have already lost their angular mo- 
mentum and formed a central cusp of compressed gas. In 
runs with higher star formation rates (e.g. the canonical 
run II26SPM), the gas clumps are partially converted to 
stars before accreting so the cusp is primarily composed 
of stars, but the net result is the same - the clumps lose 
angular momentum and form a centrally peaked rotation 
curve. 

Note that it is still possible that star formation ef- 
ficiency coupled with feedback may be important ~ in 
particular these results are consistent with the idea that 
low star formation rates combined with sufficient feed- 
back to puff up gas-dominated clumps wo uld suppress 
angula r momentum loss (indeed, the runs in Agertz et al. 
(2010) generally include feedback and always include a 
stiff equation of state, P ^ p^). However, we see that by 
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Fig. 9. — We plot the rotation curves for our canonical run (top- 
left) and the low star-formation run (top-right) at 2 = 4.5. I3oth 
exhibit a cusp of material in the core at this early epoch; however, 
in the canonical run the cusp is due to stars, whereas the low 
star-formation run's cusp is almost entirely compressed gas. This 
demonstrates that gas is first tunneled into the core of a halo even 
in the absence of significant star formation. As seen in Figure |11| 
that in the end, the stellar component ends up dominating the cusp 
by 2 = 0, regardless of how things look at 2 = 4.5. In the bottom 
two panels we also show rotation curves at 2 = 4.5 for the fixed 
physical resolution run and the 'everything' run, for comparison. 

itself, a low star formation efficiency does not change the 
distribution of matter in the simulated galaxies. 

Also shown in Figure |9] are the rotation curves at 
z = 4.5 for the fixed physical resolution run (H26SPMR), 
as well as the simulation including feedback, cooling sup- 
pression, and fixed resolution (H26SPMFCR). As can be 
seen, the cusp already appears in the fixed resolution 
run (although not quite as high as in the canonical run), 
while gas dominates in the H26SPMFCR simulation. For 
that run, the additional pressure from feedback/cooling 
suppression has allowed the infalling clumps to be incor- 
porated in the disk before losing a significant amount of 
angular momentum. For the rest of the modified physics 
runs (not shown), only the cooling suppression simula- 
tions differ significantly from the canonical run, and they 
are similar to the H26SPMFCR simulation, but with a 
somewhat more pronounced cusp. 

3.5.4. Fixed Physical Resolution 

Our canonical set of simulations used a fixed maximum 
refinement level, which translates to a best resolution 
achievable in comoving spatial units. Another common 
prescription is to fix the highest resolution as a fixed 
physical length scale (e.g. Agertz et al.||2010 |. We carry 
out a simulation of Halo 26 in which we vary the max- 
imum allowed refinement level so as to keep as close as 
possible to a physical resolution of 425 pc. Note that 
because of the factor-of-two refinement in AMR, this im- 
plies discrete resolution changes at various times during 
the simulation. The star fo rma tion history for this sim- 
ulation is shown in Figure |10[ The lower resolution at 
early times (for example, at z = 5, the spatial resolution 
is four times worse than in the canonical run) results in 
lower densities and hence a shift of the bulk of the star 
formation to later times (as in the lowered star formation 
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Fig. 10. — We plot the star formation histories for the first three of our modified runs. These runs all use identical initial conditions, those 
of canonical run H26SPM, which is also presented here. It appears that these three test runs have little effect on the bulk star formation 
history of the galaxy, although the simulation with the depressed star formation efficiency delays massive star formation for a gigayear or 
so. 
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Fig. 11. — We plot the rotation curves for the first three of our modified runs. Just as in Figures |6] & [S] we represent the circular velocity 
of the gas due to the dark-matter component (green dashed line), the stellar component (red dot-dasned line), the gas component (blue 
dotted line) and the total of all halo components (black solid line). Like Figure [Tol there is little improvement between these modified runs 
and their control run H26SPM. In fact, decreasing the star formation efficiency as in H26SPML actually drives more material into the core 
by 2 = 0. 
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efficiency run H26SPML). The rotation curve is shown in 
Figure ml and again, there is not a significant change in 
the distribution of mass in the core. 

3.5.5. Supemovae Feedback 

As described in the methodology section (Section [2]), 
we also performed one run with thermal feedback from 
Type II SN, using a moderate value of epB, as given in 
Table [l] Figure 12 shows the star formation history for 
this run, and demonstrates that the inclusion of thermal 
feedback does have some effect. There is a reduction in 
the overall burst of star formation at high redshift rela- 
tive to the canonical run with the same initial conditions. 
At late times the star formation rates are similar. How- 



ever, in Figure 13 the end result for the rotation curve 
at z = is the same as when feedback is not included. 
This occurs because the feedback is not strong enough 
to destroy the infalling clumps and prevent their loss of 
angular momentum. It is possible that stronger feedback 
will change this; however, we delay a systematic exami- 
nation of feedback prescriptions for a later paper. 

3.5.6. Cooling Suppression Models 

Many current studies employ a cooling suppression 
scheme in order to prevent cooling for a short period 
after a newly formed star is born, allowing the thermal 
feedback to efficiently operate. We performed three sim- 
ulations which integrated cooling suppression: one with 
cooling suppression alone, one with cooling suppression 
in addition to thermal feedback, and one with cooling 
suppression, thermal feedback and fixed physical reso- 
lution. Figure [l2] shows the star formation history of 
these three simulations. Interestingly, with just the lone 
addition of cooling suppression (red dot-dash line), our 
star formation is highly suppressed at early times, and 
extends to much later time, also being less bursty com- 
pared to the canonical run. The addition of thermal 
feedback on top of that further intensifies the effect of 
cooling suppression and the star formation rate becomes 
almost constant throughout the simulation at about 10- 
12 Mo/year. Finally, the synthesis of feedback, cooling 
suppression and fixed resolution results in an even lower 
overall star formation rate. In this last run, we see a step 
function in the SFR at two redshifts which correspond to 
the epochs when we change the allowed maximum refine- 
ment level, in order to preserve fixed resolution. These 
transitions show up particularly clearly in this case be- 
cause the cooling suppression operates only in the local 
cell, and since most of the disk is refined to the maxi- 
mum level, when this changes, it drastically lowers the 
efficacy of the cooling suppression, leading to more star 
formation. This demonstrates the sensitivity of our star 
formation results to the chosen parameters of the cooling 
suppression model. 

In Figure [TSj we show the resulting rotation velocity 
curves. We finally have an effect that has a significant 
impact on the rotation curves. Even cooling suppression 
by itself results in a peak rotation rate which is 100 km/s 
lower than the canonical run. Adding thermal feedback 
on top of this decreases the peak even more, to just over 
400 km/s, and adding a fixed physical resolution brings it 
down to about 350 km/s. This occurs because the cool- 
ing suppression model (and the effective feedback that it 
permits) acts to decrease the density of infalling clumps. 



These clumps are then disrupted before they can lose a 
significant amount of angular momentum and so end up 
rotating at larger radius then they otherwise would have. 
We speculate that cooling suppression operates even in 
the absence of feedback because other forms of heating, 
such as shock heating and adiabatic compression, can 
play an important role. We note that the resulting ro- 
tation curves, while much more realistic, are still some- 
what too strongly peaked at small radius, a characteristic 
shared (to a larger or smaller degree) by essentially all 
simulations that include cooling suppression. 

Finally, Figure [14] displays visual renderings of the 
gaseous and stellar components of the halos, similar to 
those of the canonical runs in Figure [4j The simulations 
which successfully reduced the buildup of material in the 
core of the galaxies have larger cold gas densities. There 
is still a lot of material in the core, but much of it is 
in compressed cold gas instead of stars. In the star pro- 
jections, young stellar disks are present and embedded 
in halos of older stars. These disks have no significant 
bars, but are coaligned with their gas disks, as we might 
predict for systems of this type. 

4. DISCUSSION 

These results demonstrate that it is challenging to gen- 
erate disk systems with the correct mass distribution. 
Without effective feedback, the default outcome is for 
dense clumps to lose angular momentum and result in 
centrally-cusped rotation curves. In particular, the sim- 
ulation with a very low star formation efficiency nicely 
demonstrates that this result is fundamentally a dynam- 
ical one, and does not depend on whether the clumps are 
primarily gas, or mostly stars. As long as they are con- 
centrated, they will lose angular momentum and hence 
rotational support. Although this result is not new, 
and there is a long history of SPH simulations which 
found this result much earlier, we show it here clearly 
and systematically using a completely different numeri- 
cal method (AMR). Therefore, the result is quite general. 

In addition, we went on to systematically vary our nu- 
merical parameters and investigate a range of resolution 
and feedback methods. We confirmed that only cooling- 
suppression feedback models are capable of significantly 
changing the mass distribution and hence the rotation 
curve. 

Cooling suppression models do effectively enhance 
feedback, although it is unclear h ow ph ysically mean- 
ingful this technique is (see Section 2.1.7), and a better 
approach might be to use hig h-resolution local models to 
generate subgrid models (see Yepes et al.|1997 Tasker &; 



|Bryan|2006[|Ceverino fc Klypin|!i0001|Joung et al.|!iOOOb 

for some attempts in this direction). 

Another constraint is the baryon content of galactic 
halos: a variety of techniques have been used to infer 
that the baryon-to-dark-matter rati o in galaxies is much 



smaller than the cos mic mean (e.g. Moster et al. 2010 



Behroozi et al.|2010 1 , implying that a significant amount 



of mass has been ejected from galactic systems (or never 
accreted in the first place). For example, Milky- Way 
massed halos only appear to host 20% of their baryons, 
with t he fraction decreasin g rapidly for smaller-mass sys- 
tems ( Behroozi et al.|20T0 ) . We find that all of our simu- 
lations result in very high disk baryon fractions; even the 
run with feedback and cooling suppression has a baryon 
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Fig. 12. — We plot the star formation histories for the last four of our modified runs. These runs all use identical initial conditions, 
those of canonical run H26SPM, which is also presented here. While feedback and cooling suppression each individually have some effect 
on lowering overall star formation, together their effects are amplified to keep star formation low throughout the galaxy's lifetime. The 
addition of a fixed physical resolution lowers that star formation even further. Unfortunately, because of the way we have implemented a 
"fixed" physical resolution in Enzo, there are discrete jumps in the star formation rate at times when we increase the comoving resolution 
(e.g. t ~ 2.5 & t ~ 7 gigayears). 
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Fig. 13. — We plot the rotation curves for the last four of our modified runs. Just as in Figures [m and Is] we represent the circular velocity 
of the gas due to the dark- matter component (green dashed line), the stellar component (re d do t-dasned line), the gas component (blue 
dotted line) and the total of all halo components (black solid line). In agreement with Figure [T2| we see that when working in concert the 
effects of feedback and cooling suppression are intensified to dampen the dense stellar cusp m the interior of our halo. The addition of 
the fixed physical resolution further decreases its rotation curve cusp and makes it look much more akin to something we might find from 
observations. 
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fraction of about 63% of the cosmic mean. This appears 
to be a general issue with cosmological galaxy simula - 
tions (see also the discussion in f Avila- Reese et al.||2011 ). 

We demonstrate that simulations must be run to z ~ 
in order to gauge the efhcacy of model parameters at min- 
imizing the effects of t he angular momentum problem. 
Many past s tudies (e.g. Ceverino fc Klypiii||2009 Joung 
et al.|[2009a ) have traded ott simulation ruri time tor in- 
creased resolution in their simulations. While there are 
some early indicators of the angular momentum problem 
like bursts of star formation and peaked rotation curves 
at redshifts as early as z = 5, successful results during 
this epoch do not guarantee successful results at z ~ 0. 
We specifically demonstrate this in the case of the low 
star formation efficiency run H26SPML, which staved off 
early bursts of star formation but eventually succumbed 
to the same fate as runs with a normal star formation 
efficiency. 

4.1. Comparison to Previous Work 

In agreement with previous work using SPH, we find 
that unless we include efficient feedback, the resulting 
systems are dominated by a too-large spheroidal compo- 



center (e.g. Navarro & Benz||1991 


iNavarro & Steinmetz 


1997 


Weil et al.||1998; Abadi et al 


tallD'ongkia etal. 


2006 


Zavala et al.||2008^ |Piontek & Steinmetz||2011|). 



TTnumber of reasons have been suggested for this in the 
past, including purely numerical causes, such as angular 
mo mentum transfer betwe en the disk and the hot halo 
(e.g. Okamoto et al. 2005[), or between cold gas c lumps 
and the hot halo (e.g. |'l'hacker fc Couchman||2000D . The 
concern was particularly that yPH simulations might be 
susceptible to this issue because of smoothing between 
hot and cold phases. However, the fact that AMR sim- 
ulations - which use a completely different numerical 
method to solve the fluid equations - find the same re- 
sult, is an indication that these effects do not dominate. 

Another suggested source of numerical angular mo- 
mentum loss is in the form of gravitational instabili- 
ties which arise from inad equate resolution of the Jeans 
lengt h in the disk (Truelove et al.|19"97{ Robertson et al. 
2004 1 . We tested this idea by running with and with 



out an additional numerical ( "Jeans" ) pressure designed 
specifically to ensure that the Jeans length was ade- 
quately resolved, finding no difference in our results^ 



A third num erica l reason is the lack of resolution ( Gov- 
ernato et al. 2004 Kaufmann et al.||2007[ ); however, we 
specifically test this over the computational range avail- 
able for us, and find no significant difference from 200- 
1 700 pc. This is in agreement with the SPH simulations 
of Piontek & Steinmetz (20 111), who also varied their nu- 
merical resolution over a wide range, and found no dif- 
ference. 

There have been a number of recent cosmologica l AMR 
simulations which we can compare to JgdiQet ar] ( |2010[ ) 



used the ART code (Kravtsov et al. 



I997p to simulate a 



halo which is smaller by an order of magnitude (about 
10^^ M0), finding peaked rotation curves (decreasing as 
the star formation density criterion was reduced). Al- 
though the halo masses are quite d ifferent, the essential 



result seems to be in agreement. Agertz et al. (20101 



used the RAMSES AMR code ('Teyssier 2002 ) tosmT 



to that found here. They argued that a low star forma- 
tion efficiency by itself was enough to produce nearly flat 
rotation curves (and that feedback was only efficient if 
extreme amounts of energy were injected). We have not 
been able to confirm the first suggestion ~ using a range 
of low efficiencies for star formation, we find that clumps 
lose angular momentum at high-redshift, and generate 
steep rotation curves, whether they are in gas or stellar 
form. The efficiency controls when the gas is converted 
to stars, but has little impact on the distribution of the 
material. In this, we are in agreement with previ o us SPH 
work (e.g.lWeil et al.|1998[ [D'onghia et al.|2006|p^iontek 
' fc Steinmetz||201l[ ) 



i"'inally, we have found that the only way to signifi- 
cantly decrease the peak of the rotation curve was to in- 
troduce a sub-grid model which enhanced the efficiency 
of stellar feedback. We briefly explored the cooling sup- 
pression model and found this to be effective. This 
agrees with a substantial nu mber of SPH simulations 



1997 



Thacker 



fc Couchman||2000 ; Sommer-Larsen et a . 2003| Stinson 


et al.||2006l |Governato et al.||2007| |Agertz et al.||2010| 


Colin et al.|2010||Piontek & Steinmetz 


2011l|Guedes etal. 


2011 


. In addition, Avila-Reese et al. 


(201l|) used AMR 


simu 


ations (with the AR'l ' code) and also found cooling 



suppression to be e ffective in obtaining appr oximately 



ulate a Milky- Way mass galaxy with similar resolution 



flat rotation curves. Ceverino fc Klypin (20091 also used 
the ART code, but with a different sub-grid model, argu- 
ing that a model in which stars are born with significant 
velocities relative to the nascent gas will feed energy into 
low-density regions, producing efficient feedback and flat 
rotation curves, although the simulation is only run to 
z - 3. 

5. SUMMARY 

In this work, we have carried out adaptive mesh refine- 
ment simulations of a sample of Milky- Way sized halos 
in order to better understand how numerical methodol- 
ogy impacts the content and structure predicted by such 
models. We selected a series of five halos ranging in mass 
from 0.5 — 1.1 x 10^^ Mq, each picked to have a quies- 
cent mass accretion history over the last 10 Gyr, with 
no major mergers, in order to focus on halos that have a 
high chance of hosting disk-like galaxies. We simulated 
the halos with high mass and spatial resolution, includ- 
ing gas, radiative cooling, star formation, and in some 
runs, feedback. We also took one halo and resimulated it 
a number of times, modifying the simulation parameters 
for each run. Our primary results are presented below. 

• We find that, without any sort of feedback, all 
five halos produce rotationally supported gas disks. 
They are all, however, dominated by a massive and 
concentrated stellar spheroid, resulting in rotation 
curves that peaks at about 500 km/s in the cen- 
tral few kpc. Therefore, we confirm previous SPH 
simulation work that also found that dissipation al- 
lowed dense clumps to lose angular momentum and 
produce halos which are too cuspy. 

• For one halo, we vary the comoving spatial reso- 
lution from 1700 to 212 pc and find the resulting 
cuspy halo to be completely robust against reso- 
lution change. Using constant physical resolution. 
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Fig. 14. — We render faceon and edgeon views for the last four of our modified halos at z = 0. The top two rows show each galaxy in 
gas, whereas the bottom two rows display each galaxy's stellar component. Each postage stamp has a width of 25 kpc. For the gas, edgeon 
views show a volume rendering with isocontours of the gas density, whereas faceon views are column density bricks. For the stars, we use 
the column density of stars in the faceon and edgeon views respectively. The color of a star particle is representative of its age where blue 
represents young stars and red represents old stars according to the color bars on the right. Notably, these modified runs produce young 
stellar disks and have a much higher density of gas in their cores than the canonical runs. 



rather than constant comoving resolution also has 
no effect in reducing the central concentration of 
mass for halos in the absence of other modified 
physics. 

Adding a model of local thermal feedback from 
Type II SN also has little effect, somewhat sup- 
pressing star formation, but with no impact on 
the z = rotation curve. Similarly, varying the 
star formation efficiency had little effect later than 
z ^ A. The use of an artificial pressure to ensure 
that the Jeans length was always resolved also re- 
sulted in little change to the star formation history 
or to the mass distribution. 

The only modification that did have a significant 
impact on the overall mass distribution in our ha- 
los was to suppress cooling in the vicinity of young 
stars (more precisely we suppress cooling in the lo- 
cal cell for 50 Myr after a star formation event). 
This, combined with the addition of thermal feed- 
back and resolving to a fixed physical resolution, 
led to a large decrease in the peak of the rotation 
curve (although still larger than observations indi- 
cate in the inner few kpc). 



To date, the most effective means for staving off the 
angular momentum problem is to employ cooling sup- 
pression in the vicinity of star formation events in order 
to intensify the effects of feedback. Most of the recent 
researchers in this field employ it by default, while at 
the same time they investigate the effects of other simu- 
lation parame ters on preventing the angular momentum 



problem (e.g. Agertz et al. 2010 Guedes et al. 2011) 



Our results indicate that the other parameters are sec- 
ondary to cooling suppression, and they do not seem to 
work effectively in the absence of it. We realize that the 
use of cooling suppression is mostly motivated by its ef- 
fectiveness in this regard, but we continue to search for 
other feedback parameterizations which are more physi- 
cally derived. 

There are several potential alternatives to cooling sup- 
pression as a form of feedback in cosmological hydrody- 
namics solutions of galaxy formation. Radiative feed- 



back shows some promising results (Kim et al. 20111, 
although scaling its effects to a large number of particle 
sources is a current computational challenge. Another 
option is that the cosmic rays produced by supernovae 
could be used as a means for transporting energy and 
momentum to the surrounding medium as an additional 
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fluid in a simulation (e.g., [Miniati 2001 Jubelgas et al. 
2008 1. Alternatively, modiiying the gas in dense re gion's 
to have a stiff equation of state ( Agertz et al.|20l"0 l may 
help "puff" up early collapsing pockets of gas. Perhaps 
there is even a redshift-de pendent feedback prescr i ption 



simi lar to those used by Sommer-Larsen et al 
andjOkamoto et al. (2005), where it was assume^ 



20031 
that 



the ImJ^' was top-heavy m the distant past resulting in 
more supernovae and a higher feedback efficiency in that 
epoch. We are currently investigating some of these op- 
tions to be presented in a forthcoming paper. 
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